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The Kramers-Kronig (KK) algorithm, useful for retrieving the phase of a spectrum based on the known spectral 
amplitude, is applied to reconstruct the impulse response of a diffusive medium. It is demonstrated by a simulation 
of a ID scattering medium with realistic parameters that its impulse response can be generated from the KK 
method with high accuracy. 

PACS: 42.25.Dd, 67.80.Mg, 42.25.Fx and 66.10.Cb 



I. INTRODUCTION 

Recently there has been a growing interest in light 
propagation in diffusive or turbid media [1]. This can 
be attributed to several areas of application-driven re- 
search. Diffusive media are ubiquitous in our environ- 
ment, and imaging through them is always a challenge. 
Clouds, mist, fog, dust and smoke decrease the visibility 
on land while surface waves and turbid water reduce vis- 
ibility at sea. In the medical field, the main obstacle that 
a physician encounters while trying to diagnose internal 
organs is the diffusivity rather than the absorption of the 
biological tissues. As a consequence, many methods were 
developed to image through and to investigate diffusive 
media [1]. 

One of the most intuitive methods to investigate a dif- 
fusive medium is to measure its' impulse response by a 
fast detector [2-6], since, in principle, the impulse re- 
sponse carries all the optical information of any linear 
medium. In practice, fast detectors measure the inten- 
sity of the impulse response. The 'first-light' component 
of the impulse response carries the ballistic information of 
the medium. This information can be used in the recon- 
struction of the ballistic image of the diffusive medium. 
Obviously, since the number of ballistic photons is ex- 
tremely small, the image can be partially reconstructed 
by the quasi-ballistic (or "snake") photons, which would 
impair the spatial resolution of the image. 

Assuming that the amount of detected quasi-ballistic 
photons is sufficient for the image reconstruction, in order 
to obtain high resolution images the detectors must be 
very fast to distinguish between the (quasi) ballistic and 
the diffusive photons. As a consequence such an imaging 
technique requires complicated and expensive equipment. 

This is one of the main motivations for developing 
spectral techniques [7-10], i.e., techniques which work 
in the spectral domain instead of the temporal one [4-6] . 
In one of these techniques the spectral response of the 



medium H (lo) — A (lo) exp [iip (lo)] is measured for each 
wavelength in a wide spectral range [11-13]. While an 
amplitude A (lo) measurement is a relatively simple task 
(and fast, since it does not require long averaging), a 
phase ip (lo) measurement is more complicated since it 
usually requires interferometric processes that are inher- 
ently complicated and susceptible to fluctuations. There- 
fore, phase measurements are of limited value when the 
medium varies in time. However, if the phase is recon- 
structed from the amplitude measurement and the pro- 
cess is sufficiently fast, the temporal impulse response 
can be easily generated even for a varying medium by a 
simple inverse Fourier transform. Since in most cases the 
diffusive medium can be regarded as linear (and hope- 
fully time-invariant) system the Kramers-Kronig (KK) 
method [14,15] can be used to derive the phase from the 
amplitude measurement. 

Nevertheless, the KK method has several drawbacks. 
Firstly, in order to derive the phase one needs the am- 
plitude over the entire spectrum (from zero to infinity), 
while experiments can produce amplitude measurements 
for only a finite range of frequencies. As a consequence, 
any evaluation of phases with the KK method is always 
an approximation (see, for example, ref. [16]). Some 
methods were developed to improve these approxima- 
tions (by improving the convergence of the integrals) such 
as the singly and multiply subtractive KK relations ( [17] 
and [18] respectively). 

A more complicated problem arises from the logarithm 
function, which diverges whenever the spectral response 
of the medium H (lo) vanishes. As a consequence the full 
KK relation is [19] 
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where the P denotes Cauchy's principal value. The first 
term on the right is the ordinary KK relation, while the 
second one is the Blaschke term, which takes account all 
of the zeros uij of H (lo). 

The problem with the zeros of H (lo) is more funda- 
mental since while the first mentioned problem can be 
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mitigated, at least in principle, by measuring the am- 
plitude for a larger spectral range, the zeros of H (w) 
cannot in general be deduced from the amplitude mea- 
surements alone. This problem is confronted mainly in 
reflection measurements, where H (u>) may vanish at cer- 
tain complex frequencies (for a discussion on this subject 
see refs. [20] and [21]). However, a diffusive medium in 
the transmission configuration is a good candidate for the 
implementation of the KK method in the following case. 

In any diffusive medium the diffusion coefficient of the 
medium D and its size L determines the smallest fre- 
quency scale Suj ~ D/L 2 , which is related to its recip- 
rocal parameter to ~ L 2 / D - the mean time a photon 
dwells in the diffusive medium . If the measured spec- 
tral range Aw = w max — ui m i n is considerably larger than 
Suj then the main features of the impulse response can 
be retrieved from the KK relations. Moreover, we will 
demonstrate that in this regime it is possible to evaluate 
the amplitude |iJ (u;) | beyond the spectral range (i.e., 
lu < Lo m - m and u> > w max ) by a certain average over its 
values inside w m i n < ui < w max . In what follows we will 
revert to wavenumbers instead of frequencies, but the 
transition between the two is trivial. 

In most previous works the desired parameter was the 
refractive index of a medium, so that the KK method 
was mostly implemented for cases where the attenua- 
tion was caused by absorption rather than by scatter- 
ing. As a result, using the KK method as a tool to mea- 
sure the impulse response of a diffusive medium was not 
common. Recently, we demonstrated [12] that the KK 
method, even in its simplest and most naive form can 
be used to reconstruct the impulse response of a Fabry- 
Perot interferometer, whose response is governed solely 
by scattering. It was shown that with relatively simple 
equipment its' impulse response can be evaluated with 
very high temporal resolution (less than 200fs). 

In this paper we demonstrate by a numerical simula- 
tion that the KK method can be useful in determining 
the phases of the transfer function of a diffusive medium. 
It is also shown that the calculated spectral response can 
be used in reconstructing the medium's impulse response 
with high accuracy. 



of scattering coefficients are similar the fine structure of 
each scatterer is not important; similarly, the exact lo- 
cations of the scatterers have a negligible effect on the 
diffusion coefficient). 

For simplicity we ignore polarization effects, and thus 
the ID stationary-state wave equation has the form 
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where ip (x) is the electromagnetic field and n (x) is the 
refractive index. 

In general, the refractive index of a diffusive medium, 
modeled as a homogenous medium with multiple ran- 
dom scatterers, has the general form presented in fig.l. 
For simplicity, however, we choose to simulate the vari- 
ations in the refractive index of the medium by delta 
functions, which corresponds to scatterers whose dimen- 
sions lj are considerably smaller than the beam's wave- 
length lj << A, consistent with our previous assumption. 
It should be stressed, that any small ID scatterer (with 
respect to the wavelength) can be replaced for any prac- 
tical reason with a delta function, the reason being that 
neither its width lj nor its strength Srij (the difference be- 
tween its refraction index and the surrounding) appears 
separately in the scattering solutions, only their product 
IjSrij does. Therefore, a small scatterer, whose width and 
strength are lj and Srij respectively can be replaced by 
a delta function whose prefactor is proportional to IjSrij 
(see below for details). 
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FIG. 1. Refractive index as a function of location for a 
diffusive medium 



II. THE MODEL 

We investigate a ID homogenous medium with N small 
scatterers, each having a different refractive index and 
width, randomly distributed in the medium. For sim- 
plicity it is assumed that the width of each scatterer lj is 
considerably smaller than the beam's wavelength A (i.e., 
lj « A for every j), however, this is not a restrictive as- 
sumption since a diffusive medium, whose dimensions are 
considerably larger than the medium's diffusion length 
is characterized, almost by definition, only by a median 
value for the scattering coefficient (i.e., as long as the set 



The square of the refraction index can be separated 
into homogenous (tIq) and varying (2rioAn) parts 
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With the definition of the wavenumber k = 2nn a /\ 
the wave equation can be written 
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where 2An (x) /no — a j$ ( x ~ A?)' N is the number 

3 

of scatterers, Lj = Lj_\ +a 3 - = a m is the position of 

the jth scatterer, i.e., a 3 - are the distances between two 
adjacent scatterers, and 
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is the strength of the jth scatterer, where Snj and lj are 
the change in its refractive index and its size respectively. 

Therefore, the wave equation for the diffusive medium 
reads 
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FIG. 2. yl schematic presentation of the diffusive medium. 
Each scatterer is presented in the model as a delta function 
change in the refractive index. 

In the appendix we elaborate on the derivation of the 
medium's spectral response H (fc) from eq.6. We then 
apply the KK method to a finite spectrum fc m i n < k < 
fcmax to determine the phase: 
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which is then substituted into: 



Hkk (fc) = |ff (fc)|exp [*¥>**(*;)] 



(7) 
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to evaluate the medium's spectral response. If we keep 
Afc = fcmax — fcmin >> Dn/ (cL 2 ) (where n is the refrac- 
tive index of the medium) then most of the features of 
the impulse response can be derived. Moreover, in this 
regime, approximation (7) can be improved by noticing 
that the mean value of In \H (k)\ does not change much 
beyond the measured region fc m j n < k < fc max (at least 
in its spectral vicinity). Therefore, the boundaries of the 



integral (7) can be taken as and oo while In \H (k)\ for 
< k < fc max and fc max < k < oo is replaced with its 
average value, i.e., 
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where (In \H (k)\) = (fc max - fc^)" 1 dk' In \H (k')\ 

is the mean value of In \H (k)\ in the measured region. 

It should be noted that in cases where the KK method 
is used to derive the refractive index (as is the case in 
negligibly scattering media) the variations in \n\H (k)\ 
are of the same scale as the spectral range Ak, i.e., 
d\n\H (k)\ /dk ~ Afc -1 . Therefore, in these cases an- 
alytical continuation and extrapolations are used to ap- 
proximate In | if (fc)| beyond the measured region [22]. In 
the scattering medium case the situation is considerably 
different, namely d In | H (fc)| /dk » Afc -1 , and there- 
fore extrapolations has little value. On the other hand 
the spectral variations are so strong that they rapidly 
converge to the average value. 

By solving the integrals one obtains an evaluation of 
the phases from the amplitude measurements 
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is a correction term. 

This phase is then substituted into eq. (8) to recon- 
struct the full transfer function. 

If the spectrum of the input signal is a rectangu- 
lar function in the spectral domain fc m j n < k < fc max 
(and therefore the field is proportional to _Ej„ (t) oc 
exp (ikct) sine (Afcct/2), where Afc = fc max — fc m ; n and 
k = (fc max + fc m in) /2) the exact impulse response, which 
will be measured at the end of the medium, i.e., at 
x > Ljy, is 



I(t) 



dk'H (fc') exp (ik'ct) dk' 



(12) 



while the KK reconstruction, which is based only on the 
amplitude (or intensity) measurements, predicts 
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We will demonstrate by a simulation that in a realistic 
case the latter expression is a very good approximation 
to the former. 



III. SIMULATION 

For simplicity we choose air as the homogenous 
medium, i.e., uq = 1 (however, the results can easily 
be scaled to any dielectric medium with an arbitrary n ) 
with 150 small scattcrers, where the distance between 
them is a random variable, distributed uniformly between 
and 5/xm (i.e., < aj < 5^im). Similarly, the strength 
of each scatterer is also a uniform random number so 
that < aj < 0.03^to (this corresponds to glass parti- 
cles having widths on the order of tens of nanometers). 
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FIG. 3. The impulse response of the medium. The upper 
plot is the exact reconstruction, while the lower one is the KK 
reconstruction. 

We assume that the incoming pulse that penetrates the 
medium has a rectangular spectral shape, i.e., E in (i) oc 
exp (ikct) sine (Akct/2), with u> = kc = 2ir x 750THz and 
Auj = Akc = 2irx 300THz. 

Since D = clR\y/ n where Irw — o,/R (a is the aver- 
age distance between scatterers and R is the mean re- 
flection coefficient) is the random walk in the diffusion 
process, then the minimum spectral range required is 
5k = a/RL 2 = (NRLy 1 (where N is the total num- 
ber of scatterers). 

In this case the measured spectral range is six orders 
of magnitude larger than 8k i.e., 5k = 15m _1 << 6.3 x 
10 6 ™- 1 Afc. 

In Fig. 3 the reconstruction of the impulse responses 
of the exact solution / (t) and the KK reconstruction 
Irk (t) are presented. 
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FIG. 4. Zoom-m of fig. 3 

As can be seen, except for the delay between the two 
pulses, which is a direct consequence of the KK tech- 
nique, the two signals are very similar. Evidently, they 
are not identical but the differences between them are 
quite small. In Fig. 4 we reveal details of the temporal 
response and compare the two responses over a small 
temporal window, showing that the two are alike even 
on the microscopic level. 
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FIG. 5. Same as Fig. 3 but with a spectral width of 
Alu = 2tt x SOTHz. 

If the spectrum of the incoming pulse is narrower, i.e., 
the incoming pulse is temporally broader, the spikes in 
the outgoing pulse are respectively wider. In Fig. 5 
and 6 the impulse response of the same system is pre- 
sented when the spectral width of the incoming pulse 
is ten times narrower (than in Figs. 3 and 4), that is, 
Alu = 2n x 30THz. As can be seen from these plots, 
the KK approximation is even better for the spectrally 
narrow pulse. Therefore, although the integration in the 
KK expressions covers a narrower region, and in principle 
should yield inferior results, since the spikes are tempo- 
rally wider and therefore are less susceptible to dispersion 
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deformation, the reconstruction is improved. The prob- 
lem is, however, that as Afc decreases and gets closer 
to 8k there is insufficient spectral information to recon- 
struct the complete impulse response and the error in the 
evaluation of macroscopic averages (such as the diffusion 
coefficient) increases. 




FIG. 6. Same as Fig. 4 but with a spectral width of 
Alo = 2tt x SOTHz. 



IV. THE EFFECT OF THE CORRECTION TERM 

To illustrate the importance of the correction term 
AifKK (k) in cq.8 we repeat the simulation for a band- 
width Aui = Akc = 2ir x 150THz with and without the 
CT. 

In the upper panel of Fig. 7 the reconstruction was 
made with the CT, and in the lower panel the CT was 
absent. There is a clear improvement in the upper panel. 




FIG. 7. In upper figure the reconstruction was done with 
the CT, while in lower one the calculations were done without 
it. In both cases the solid line represents the direct calculation, 
and the dashed line corresponds to the KK reconstruction. 

To understand the influence of the correction term we 



can expand it with respect to the zero-correction point 

^0 = y/k-miu ^max 
AfKK (k) = 

-± (In \H (k)\) [4k - 2riK 2 + f (l + §r, 2 ) k 3 + O (k 4 )] 

(14) 

where r\ = ^ is a relatively small parameter, which char- 
acterizes the normalized spectral width (Afc = fc max — 
fc m i n ) and k = ^j?- is the normalized wavenumber. 

The first term in the expansion is responsible for a 
constant time delay, which therefore has a trivial influ- 
ence on the solution. Moreover, this linear delay, unless 
\H (fc)| is very small due to the multiple scattering, is on 
the same scale as the initial pulse width r = (cAfc) -1 , 

2 
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When Afc << fc , i.e., r\ « 1 the second term, which 
is the first dispersion term, can be neglected. Therefore, 
the first term which causes the main deformation is of 
the third order, so that 



A(fi non-linear (fc) - - A (In \H (fc)|) K \ (16) 

The dispersion coefficient is therefore proportional to 
(\n\H{k)\). 

Since |k| < 1 the time-delay, which is a consequence of 
this term is bounded by 
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(17) 

which is independent of the spectral width Afc. From 
this respect, the deformation in the signal is proportional 
to the peaks' temporal width, however, when Afc de- 
creases the average (In \ H (fc)|) is a better approximation 
to the value of ln\H (fc)| outside the measured spectral 
width Afc, and therefore the KK reconstruction is im- 
proved. 



V. SUMMARY 

We have demonstrated through numerical simulations 
that the KK method can be used to reconstruct the im- 
pulse response of a scattering medium with realistic pa- 
rameters. It was shown that when the measured spectral 
width is considerably larger than the reciprocal of the 
diffusion length, i.e. Afc = fc max — fc m ; n >> Dnj (cL 2 ), 
the KK method yields a satisfactory prediction of the im- 
pulse response. Moreover, it was demonstrated that it is 
possible to take advantage of the fact that in a diffusive 
medium the spectral variations are very large but their 
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running average has relatively small variations. There- 
fore, the integrand of the KK relations can be evaluated 
even outside the measured spectral domain as the aver- 
age value of the measurements. It was shown that this 
evaluation improves the reconstruction of the impulse re- 
sponse. 

Although we focused in this work on the optical prop- 
erties of a diffusive medium, this technique can be im- 
plemented to any wave dynamics in diffusive media, e.g., 
acoustic scattering (photon scattering), x-ray scattering, 
electron scattering etc. 

Owing to the simplicity and speed of the KK method, 
we believe that this technique is a promising tool for char- 
acterizing and imaging through diffusive media. 

This research was supported by the ISRAEL SCI- 
ENCE FOUNDATION (grant no. 144/03-11.6). 
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Similarly, the discontinuity of the field at x = 
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Therefore, the influence of each of the scatterers can 
be describe by the 2x2 matrix, 



which relates Vj 



to v 



by 
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= A 3 V j 
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VI. APPENDIX A: CALCULATIONS OF THE 
SPECTRAL RESPONSE (TRANSFER) 
FUNCTION H (K) 



u + j_ { exp(ikx) u + j exp(ikx) 
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FIG. 8. An illustration of the incoming and outgoing waves 
at the vicinity of a single scatterer. 

Since we are discussing a ID medium, consisting of 
the scatterers in an otherwise totally homogenous envi- 
ronment, then in every point in space between two scat- 
terers the field of the incoming and outgoing waves (see 
Fig. 8) can be described by two coefficients. Therefore, 
for a given k the vector 



(18) 



fully describes the field between the jth and (j+l)th scat- 
terers, which can be written as 



1 — irjj irjj 
-irjj 1 + irjj 



(25) 



and rjj = ajk/2. The homogenous medium between the 
jth and (j+l)th scatterers can be described by the ma- 
trix 



e~xp(ikaj) 

exp (— ika,j) 



(26) 



With these two types of matrices, one can generate the 
matrix of the medium that includes the first j scatterers 
Mj if the matrix of the j-1 scatterers Mj_i is known 



Mj = AjDjMj-i. 



(27) 



Therefore, with a simple recursion the matrix of a 
medium Mjy with an arbitrary number of scatteres N 
can easily be generated. 

Finally, the transfer function H (fc), which in our case 
is the transmission coefficient of the medium, can easily 
be reconstructed from the matrix of the entire medium 
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by 



for every k. 



H (k) = det (Mjv) /m 22 



(28) 



(29) 



V> (xj < x < Xj + i ) = u~j exp (ikx) + u- exp {—ikx) 



(19) 



By applying the continuity of the field at the two ends 
of the scatterer at x = 



<4>(x = -0) =ip(x = +0) . 



(20) 
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